Quasi-harmonic approximation of thermodynamic properties of ice Ih, II, and III 
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Several thermodynamic properties of ice Ih, II, and III are studied by a quasi-harmonic ap- 
proximation and compared to results of quantum path integral and classical simulations. This 
approximation allows to obtain thermodynamic information at a fraction of the computational cost 
of standard simulation methods, and at the same time permits studying quantum effects related to 
zero point vibrations of the atoms. Specifically we have studied the crystal volume, bulk modulus, 
kinetic energy, enthalpy and heat capacity of the three ice phases as a function of temperature and 
pressure. The flexible q-TIP4P/F model of water was employed for this study, although the results 
concerning the capability of the quasi-harmonic approximation are expected to be valid indepen- 
dently of the employed water model. The quasi-harmonic approximation reproduces with reasonable 
accuracy the results of quantum and classical simulations showing an improved agreement at low 
temperatures (T< 100 K). This agreement does not deteriorate as a function of pressure as long as 
it is not too close to the limit of mechanical stability of the ice phases. 

PACS numbers: 65.40.-b, 65.40.De, 63.70.+h, 62.50.-p 
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I. INTRODUCTION 

Computer simulations of water and ice remain a very 
active research field after the first pioneering works more 
than 40 years agoJ^ Despite the impressive progress 
achieved in the development of total energy computa- 
tional models, simulation methods, and computer archi- 
tectures, we are still far from having a standard model 
to simulate water at the atomistic level. In fact, water 
molecules are currently simulated with various degrees 
of sophistication. In the order of increasing complex- 
ity, one finds water models as point particles with short- 
range anisotropic interactions^ as rigid molecules^ as 
flexible molecules with either harmonic^ or anharmonic 
OH bonds, 7 or as flexible polarizable molecules Accu- 
rate characterization of the electronic structure of the 
molecule and its chemical reactivity is only possible 
when ab initio methods, like density functional theory 
(DFT)^2ri& are used. Within this field, the develop- 
ment and application of new functionals specially de- 
signed to treat van der Waals interactions is a focus of 
recent interest J^r— 

An additional aspect in the simulation of water phases 
is performing the calculation of thermodynamic proper- 
ties either in a classical or a quantum limit. Evidence of 
the relevance of quantum effects related to nuclear mo- 
tions in water phases is provided by measurements of 
isotope effects, that would vanish in the classical limit. 
As an example the melting point of ice Ih under normal 
conditions is shifted by 3.8 K in deuterated ice, and by 
4.5 K in tritiated ice. Interestingly, the density of ice Ih 
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shows an inverse isotope effect, i.e., D2O ice is expanded 
with respect to H2O ice ; 20 ' 21 whose microscopic origin 
has been recently revealed^ 

The path integral (PI) formulation of statistical me- 
chanics is the most applied method to treat quantum 
nuclear motion in water] 7 ' 8 ' n i 14 i 23 ~ — However, the com- 
putational cost of this method is much larger than that of 
its classical counterpart, and for this reason many equi- 
librium properties of water and ice have not been yet 
studied by quantum simulations. A prominent example 
is the complex phase diagram of ice, with fifteen differ- 
ent phases. To date, the simulation of relevant parts of 
this phase diagram has been only accessible to classical 
simulations using effective rigid water modelsj 34 i 35 The 
computational overhead of the quantum PI simulations 
is caused by the need to work with multiple replicas of the 
system. Several approximations are available that help 
to reduce it. For the case of empirical point charge mod- 
els of water, the ring polymer contraction scheme allows 
for a significant increase of computational efficiency^ In 
the case of ab initio methods the number of replicas can 
be reduced dramatically by means of an appropriate gen- 
eralized Langevin equation^ 

The quasi-harmonic approximation (QHA) has been 
occasionally applied to study thermodynamic properties 
of ice phases. Here the collective vibrations are described 
by a set of harmonic oscillators and thus the partition 
function can be expressed as an analytic function of 
the crystal volume and the temperature. An advantage 
of this approximation is that all equilibrium thermody- 
namic properties can be derived straightforwardly. This 
implies that the computational effort of the QHA is only 
a very small fraction of that required for a PI simulation 
of an ice phase, and that this fraction is independent of 
the chosen model potential. Further advantages of the 
QHA are the absence of statistical errors, as opposed to 
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any classical or quantum simulation, and the possibility 
to account for finite size effects by a Brillouin zone inte- 
gration of the phonon dispersion curves, rather than by 
increasing the size of the cell. 

The QHA in combination with ab initio DFT has al- 
lowed the explanation of the inverse isotope effect in the 
crystal volume of ice Ih at atmospheric pressure. 22 Also 
the negative thermal expansion of ice Ih at low temper- 
atures has been studied by the QHA^ as well as the 
elastic moduli and mechanical stability of the hydrogen 
ordered ice VIIIi 3 -^ In addition, the mechanical stabil- 
ity of ice Ih under pressure has been studied by this 
approximation.— One conclusion of these results is that 
the QHA seems to be reasonably realistic for ice phases. 
However, a systematic study of its potential as a predic- 
tive tool is still missing. 

The purpose of this work is to check the ability of the 
QHA to predict thermodynamic properties of ice phases. 
We have analyzed three different phases (ices Ih, II, and 
III) in order to have a broad field of comparison. Our 
aim here is to compare PI results to the expectations of 
a much simpler QHA. We have used the point charge, 
flexible q-TIP4P/F model to investigate several thermo- 
dynamic properties as a function of both pressure and 
temperature* 7 This model has been recently employed 
in PI investigations of isotope effects) 30 ' 31 proton kinetic 
energies, 33 and pressure effects in water and ice Ih^2 We 
expect that the validity of the QHA will be to a large 
extent independent of the employed potential model. 

The layout of the manuscript is as following. Compu- 
tational details of the numerical simulations and QHA 
of the ice phases Ih, II, and III are given in Sec. [TTJ 
The phonon density of states (DOS) and the Griineisen 
constants calculated with the q-TIP4P/F model are pre- 
sented in Sec. IIII1 Comparison of QHA and simulation 
results for various pressures (P) and temperatures (T) 
are presented in Sec. IIV1 Specifically, we study the vol- 
ume (V), bulk modulus (£?), kinetic energy (K), enthalpy 
(H), and heat capacity (C p ). Also the comparison to ex- 
perimental data is presented whenever available. The 
paper closes with the conclusions. 



II. COMPUTATIONAL PROCEDURES 

In this Section we summarize some relevant computa- 
tional details concerning the ice cells, the PI simulations 
and the QHA. 



A. Simulation cells 

Ice Ih displays proton disorder, i.e., while the oxygen 
atoms occupy fixed lattice positions, the water molecules 
display random orientations with the constraint that the 
resulting H-bond network must be compatible with the 
Bernal- Fowler ice rules, 41 These rules imply that each 
oxygen is tetrahedrally coordinated to two H atoms by 



strong OH covalent bonds (~490 kJ / mol) and to two ad- 
ditional H atoms by weaker H-bonds (^25 kJ/mol). Ice 

II is an ordered phase where the molecular orientation is 
fully determined by the crystal symmetry. Ice III displays 
partial proton disorder, i.e., the occupancy of 2/3 of the 
crystallographic H sites deviates from 50% and it is found 
between 35 and 65%4 2 - Proton disordered structures with 
nearly zero dipole moment were generated for ice Ih and 

III by a simple MC algorithm designed to explore molec- 
ular orientations that obey the ice rules^ 3 - In line with re- 
cent computer simulations a full proton disordered struc- 
ture was assumed for ice Ill^dids as this simplification 
seems to have only a minor effect in the phase diagram of 
ice^ The largest ice cells employed in our study contain 
N = 288 molecules for ice Ih, and 324 molecules for ices 

II and III. The ice Ih cell was orthorhombic with param- 
eters (4ai, 3-\/3ai, 3a3), with (ai,a3) being the standard 
hexagonal lattice vectors of ice Ih^ while ice II and ice 

III were studied by 3 x 3 x 3 supercells of the crystallo- 
graphic cell, which belong to the rhombohedral and the 
tetragonal crystal systems, respectivel y 42 ' 47 Simulations 
were conducted also for smaller simulation cells to check 
the convergence of the results with the system size. In 
the following the computational conditions employed for 
the PIMD simulations are presented. 



B. Path integral simulations 

In the PI formulation of statistical mechanics the quan- 
tum partition function is calculated through a discretiza- 
tion of the integral representing the density matrix. This 
discretization leads to a suggestive picture: with respect 
to its equilibrium thermodynamic properties, the quan- 
tum system results to be isomorphic to a classical one. 
However, this isomorphism does not apply to the dy- 
namic (time dependent) properties of the system. Specif- 
ically, the classical isomorph is constructed by a simple 
substitution of each quantum particle (here, atomic H 
and O nucleus in the water molecule) by a ring polymer 
of L classical particles (beads), connected by harmonic 
springs with a temperature- and mass-dependent force 
constant. Details of this practical simulation method can 
be found elsewhere4£~— Equilibrium properties in the 
classical isomorph can be derived by a classical molec- 
ular dynamics (MD) algorithm, that has the advantage 
against a Monte Carlo method of being more easily par- 
allelizable, an important fact for efficient use of mod- 
ern computer architectures. Effective reversible integra- 
tor algorithms to perform PIMD simulations have been 
described in Refs. I52T4551 Ref. [H introduces useful 
methods to treat full cell fluctuations and multiple time 
step integration. All simulations were done using origi- 
nally developed software and parallelization was imple- 
mented by the MPI library^ The potential energy of 
each replica of the system is calculated in parallel by a dif- 
ferent processor. Other interesting optimization schemes, 
not used in this work, are based on ring polymer contrac- 
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Figure 1: Results for the convergence of the enthalpy (a) and 
volume (b) of ice II with the number of beads (L) employed in 
the PIMD simulations. The constant shifts that correct the 
finite L results derived with the condition LT = 6000 K are 
explicitly shown. The simulations were performed at T —200 
K and P =1 atm. 



tion techniques^ 

PIMD simulations in the isothermal-isobaric NPT en- 
semble (TV being the number of water molecules) were 
conducted for ice Ih, II, and III at temperatures up to 
300 K and pressures between -1 and 2 GPa. We have em- 
ployed the point-charge, flexible q-TIP4P/F model, that 
was originally parameterized to provide the correct liquid 
structure, diffusion coefficient, and infrared absorption 
frequencies in quantum simulations^ Periodic boundary 
conditions were applied to the simulation cell and the 
Ewald method was employed to calculate the long range 
electrostatic potential and atomic forces. Expected aver- 
ages were derived in runs of 10 6 MD steps (MDS) using 
a time step of At = 0.2 fs. The system equilibration was 
conducted in runs of 10 5 MDS. To have a nearly con- 
stant precision as a function of temperature, the number 
of beads L was set as the integer number closest to fulfill 
the relation LT = 6000 K, i.e., at 200 K the number of 
beads was L = 30. The classical limit is easily achieved 
by setting L—l in the PIMD algorithm. The largest er- 
ror associated to the finite number of beads is caused 
by the highest energy vibrations, i.e., the intramolecular 
OH stretching modes. These modes remain essentially 
in their ground state as their energy quantum (Hujoh) 
is several times larger (~ 16) than the thermal energy 
(fcflT) at the highest studied temperature. Thus we ex- 
pect this error to be a nearly T (and P) independent 
shift. The shift has been estimated for the quantities 
of interest (volume, enthalpy, and kinetic energy) by a 
series of NPT simulations of the three ice phases using 
different numbers of beads at T= 200 K and P= 1 atm. 



The L — > oo limit was extrapolated from the best linear 
fit in either the variable 1/L 2 or 1/L. In Fig. [T]the result 
of this extrapolation is represented for both the enthalpy 
and volume of ice II. The constant shifts that correct the 
thermal averages derived with the relation LT — 6000 
K are 3.71 kJ/mol for the enthalpy and 1.89 kJ/mol for 
the kinetic energy. The corresponding correction shift for 
the volume amounts to -0.03 A 3 /molecule for ice Ih and 

III, and -0.06 A /molecule for ice II. Additional compu- 
tational conditions in the simulations were identical to 
those reported in Ref. [30l . 



C. Quasi- harmonic approximation 

1. Basic equations 

The quasiharmonic approximation employed here is 
based on the following three assumptions: 

a) the simulation cell, defined by the vectors 
(ai,a2,a3), is considered to scale isotropically with the 
crystal volume, V, as 



a,(V) = (V/V ref y/h 



4,ref 



1,2,3,. 



(1) 



V re f is the volume of the reference cell 
(ai, r e/, a2,re/> a 3,re/) that minimizes the potential 
energy at a certain chosen pressure, P re f, where the ice 
phase is mechanically stable. 

b) the lattice vibrations are described by harmonic os- 
cillators of wavenumber Wk, with k combining the phonon 
branch index and the wave vector within the Brillouin 
zone. 

c) the wavenumbers Wfc depend on the volume of the 
crystal, which may change with temperature and pres- 
sure. However, for a given volume the wavenumbers ujk 
are constants (independent of T and P) as in a harmonic 
approximation. 

With these assumptions, the Helmholtz free energy of 
the ice cell with N molecules in a volume V and at tem- 
perature T is given by 



F(V, T) = U S {V) + F„(V, T) - TS H 



(2) 



where Us(V) is the static zero-temperature classical en- 
ergy, i.e., the minimum of the potential energy when the 
volume of the ice cell is V. The entropy Sh is related 
to the disorder of hydrogen, it vanishes for ice ordered 
phases as ice II. Sh was estimated by Pauling for fully 
disorder phases using simple counting schemes of allowed 
water orientations^ 
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F V (V,T) is the vibrational contribution to F, 



F V (V,T) = Y1 



P 



In [1 - exp (-fihujk)] 



(3) 



(4) 
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Here /? is the inverse temperature: l/fc^T. This expres- 
sion for F v is valid for quantum harmonic oscillators. If 
one is interested in the classical limit of the QHA the 
vibrational contribution changes to 



(5) 



The terms Us{V) and Sh remain the same for either a 
quantum or classical limit. The calculation of those ther- 
modynamic properties derived from the partition func- 
tion (as the Gibbs free energy, enthalpy, internal energy, 
kinetic energy, state equation, thermal expansion, bulk 
modulus, heat capacity, etc.) can be easily performed 
with the help of Eq. ©. For example, the Gibbs free en- 
ergy, G(T, P), can be derived by seeking for the volume, 
Vmin, that minimizes the function F(V, T) + PV, as 



G(P,T) = F(V min ,T) + PVrr, 



2. Direct implementation 
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The direct implementation of the QHA implies the fol- 
lowing steps: 

a) determination of the ice reference cell by an energy 
minimization where both water molecules and cell pa- 
rameters are simultaneously optimized. Both the shape 
of the reference cell and its volume V re f are obtained in 
this step. 

b) a set of cell volumes {V^} of interest is selected. 

c) given a cell volume V,, the simulation cell is derived 
by the scaling of the reference cell as shown by Eq. ([1} . 
The water molecules are then relaxed to their minimum 
potential energy and the crystal phonons are calculated 
and tabulated. This tabulation allows us the calculation 
of F(Vi,T) for the set of volumes {V,} by means of Eq 
©■ 

d) for the purpose of finding the minimum of F(V,, T) 
as a function of V (or alternatively, the minimum of G = 
F + PV), it is convenient to perform a polynomial fit 
of F(Vi,T) and then look for the minimum of the fitted 
function. 

This implementation requires, for each of the studied 
volumes Vi , a crystal structure minimization and the cal- 
culation of the corresponding crystal phonons. Typically 
we have selected a set of 50 different V, for each ice phase 
and chosen a 5th degree polynomial as fitting function to 
find the minimum of F(V, T) as a function of V. Using 
an empirical potential this implementation is straightfor- 
ward and does require a little amount of computer time 
even for simulation cells containing several hundreds of 
water molecules. 

The crystal phonon calculation has been performed by 
the small-displacement, or frozen phonon methodj 58 i 59 
The basic principle is to compute the force constants 
between atom pairs numerically deriving the (analytic) 
forces. Atoms are displaced a small, but finite, amount 
from their perfect-lattice positions, and all the atomic 



forces generated by this displacement are calculated. 
Given that the displacement is small, the force constant 
describes the proportionality between displacements and 
forces. The atomic displacement employed in this work 
is Sx — 10~ 6 A along each cartesian direction. Although 
we have used a flexible water model, the method can be 
also applied to rigid water models by allowing molecu- 
lar displacements associated with small translations and 
rotations of the rigid units. See Ref. [6(| for a full ac- 
count of the calculation of the external phonon modes 
associated to rigid molecules. 

The QHA results have been derived by a V (k = 0) 
sampling of the Brillouin zone of the ice simulation cell. 
Although this is an adequate approximation, given the 
large size of the employed cells, its main justification is 
that we are interested in the comparison to PIMD simu- 
lations. Note that the application of periodic boundary 
conditions in the simulation implies that all atomic im- 
ages move in phase, and this corresponds exactly to the 
spatial symmetry of the k = phonons. 



3. A simpler implementation 

For calculations based on non-empirical potentials of 
water, e.g., in ab initio DFT theory, it may be computa- 
tionally convenient to reduce the number of energy mini- 
mizations and phonon calculations by using the phonon- 
dependent Griineisen parameters defined as^ 



7fe 



d\nV 



(7) 



Vref 



This equation can be integrated assuming that jk is a 
volume independent constant to obtain the volume de- 
pendence of Wfe as 



Wfc(F) « uj k (V re f) x 



V 



V 



ref 



-Ik 



(8) 



A linear expansion of the previous relation, i.e., taking 
the derivative (duik/dV) as a constant independent of V, 
leads to the alternative expression 



UJ k (V) « U k (V r ef) X ^ 



1 



v-v 



ref 



V 



ref 



Ik 



(9) 



The last two relations allow us to calculate the func- 
tion F V (V,T) in Eq. Q from a tabulation of the values 
(wfc,7fc) obtained for the reference cell. We have com- 
pared the direct implementation of the QHA (see Sect. 
Ill C 21) with the simpler implementation of Eqs. (|SJ) and 
© by calculating the function V(T) at P = 0. We find 
that Eq. © works slightly better than Eq. © for ice 
Ih, while for ice II and III the opposite behavior was 
observed. A recent QHA study of the thermal expan- 
sion of ice Ih using DFT was based on Eq. ©^ Note 
that the numerical calculation of the Griineisen constants 
requires at least two energy minimizations and phonon 
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Table I: Volume, potential energy, and simulation cell param- 
eters of the reference cells of ice Ih, II, and III obtained by 
energy minimizations with the q-TIP4P/F model at P re f = 0. 
The direct implementation of the QHA implies a phonon cal- 
culation at 50 different volumes uniformly distributed in the 
interval [V min ,V m ax]- 



0.02 - — 





Ih 


II 


III 


N (molecules/cell) 


288 


324 


324 


Vref (A 3 /molec.) 


30.96 


24.14 


24.99 


U s ,ref (kJ/mol) 


-61.98 


-60.84 


-60.86 


AU s ,ref (kJ/mol) 


-1.15 





-0.02 


Ol.re/ (A) 


17.78 


22.98 


19.68 


«2,re/ (A) 


23.09 


22.98 


19.77 


a 3 ,ref (A) 


21.76 


22.98 


20.80 


a n 


90.0 


113.2 


89.9 


n°) 


90.0 


113.2 


89.9 


7(°) 


90.0 


113.2 


89.8 


V m i„(A 3 /molec.) 


29.47 


21.75 


22.48 


V m ax (A 3 /molec.) 


35.05 


27.31 


28.22 



calculations (at V re f and at an additional volume in the 
proximity of V re f). Therefore the computational cost of 
this simpler implementation is reduced by about a factor 
of 25 in comparison to the previous direct implementation 
of the QHA. In the next Section the direct implementa- 
tion of the QHA is applied to the study of ice phases. 



III. Q-TIP4P/F ICE PHONONS 

Reference cells for the QHA have been derived at 
Pref=0 for the three ice phases by energy minimiza- 
tions implying optimization of both atomic positions and 
simulation cell. Optimized simulation cells as well as 
the corresponding volumes (V re f) and potential energies 
(Us,ref) are summarized in Tab. |U The direct imple- 
mentation of the QHA involved the calculation of the 
vibrational modes for 50 different volumes uniformly dis- 
tributed in the interval \V m in,V ma x] (see Tab. |T|. Note 
that the largest volume per molecule corresponds to ice 
Ih (28% and 24 % larger than those of ices II and III, 
respectively). By setting the potential energy of ice II 
as the zero of the energy scale, the structure of ice Ih 
results more stable by A£7s jre /=-1.15 kJ/mol, a value in 
good agreement to that of -1.17 kJ/mol reported in Ref. 
1441 However, for ice III we find AUs,ref = -0.02 kJ/mol, 
significantly different from the previously reported result 
of -0.1 kJ/molJi We have checked that this difference 
is a consequence of the disorder of hydrogen in ice III. 
Support for this is given by the equilibrium potential en- 
ergies Us, ref of a series of five randomly generated ice III 
structures. These energies scatter in an energy window of 
0.2 kJ/mol, a value more than two times larger than the 
energy difference between our results and the one Ref. 
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Figure 2: Vibrational density of states obtained for the ref- 
erence cells of ice Ih, II, and III with the q-TIP4P/F water 
model. The harmonic modes correspond to a T sampling of 
the Brillouin zone and each mode was plotted as a normalized 
Gaussian with a width of 10 cm" 1 . 



44,. Interestingly this effect is much smaller in disordered 
ice Ih, where the energy fluctuations computed for five 
randomly generated structures are reduced by more than 
an order of magnitude to about 0.01 kJ/mol. 

The vibrational DOS calculated with the q-TIP4P/F 
model for the reference cells of the three ice phases are 
shown in Fig(21 The wavenumbers of the vibrational 
modes are separated into four groups comprising 3iV 
translational, 3N librational, N bending, and 2N stretch- 
ing modes. Translational and librational modes are re- 
lated to the intermolecular interactions and can be dif- 
ferentiated not only by their respective frequency win- 
dows but also by their effective masses. The translational 
modes display a "heavy" molecular mass while the libra- 
tional modes display a "light" hydrogen mass, a fact that 
can be experimentally observed by the frequency shifts 
found after isotopic substitution of either O or H atoms. 
Within a harmonic approximation this shift scales as the 
squared root of the effective mass. The intramolecular 
modes (bending and stretchings) display also a "light" 
hydrogen mass. It is interesting to observe the anticor- 
relation between librational and stretching modes in the 
ice phases: the lower the frequency of librational modes, 
the higher the frequency of the stretchings modeSi 62 i 63 
The calculated Griineisen constant for ice Ih and III are 
displayed in Fig. [3] 

A fingerprint of the vibrational structure of the ice 
phases is presented in Tab. HH that summarizes the re- 
sult of averaging the complete set of 9N wavenumbers 
calculated for the ice reference cell into 9 groups of N 
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Figure 3: Griineisen constants of ice Ih and III calculated for 
the q-TIP4P/F model. Results for ice II lie in-between ice Ih 
and III and are not shown for clarity. 



modes taken from a list ordered by increasing wavenum- 
ber. The resulting average frequencies and Griineisen 
constants reveal significant differences between the ice 
phases. The negative Griineisen constant associated to 
the 7Y translational modes of lowest energy is the origin 
of the negative thermal expansion experimentally found 
in ice Ih at low temperature.— For ice II the averaged 
Griineisen constant of these modes is also negative but 
much smaller in absolute value, while for ice III is pos- 
itive. Negative Griineisen parameters are also found for 
the 2N stretching vibrations with absolute values in the 
order Ih>II>III. We note that the Griineisen constants of 
the stretching vibrations in the q-TIP4P/F model are a 
factor of two smaller than those obtained in the ab initio 
DFT study of ice Ih. 22 These larger DFT absolute values 
are crucial to correctly describe the experimental inverse 
isotope effect found in the crystal volume of normal and 
deuterated ice Ihj22 an effect that is not reproduced by 
the q-TIP4P/F models 



IV. QHA RESULTS 

In this Section we evaluate the performance of the 
QHA in the description of ice Ih, II, and III by compar- 
ing it to PIMD simulation results. Our study is focused 
on the temperature and pressure dependence of several 
ice properties as crystal volume, bulk modulus, kinetic 
energy, enthalpy, and heat capacity. Both quantum and 
classical limits have been studied, and compared to ex- 
perimental results whenever they are available. 



Figure 4: Thermal expansion of ice Ih, II, and III at atmo- 
spheric pressure calculated with the q-TIP4P/F water model. 
The QHA results (lines) are compared to PIMD and MD sim- 
ulations (symbols). Both quantum (q) and classical (cla) lim- 
its are shown. Error bars of the simulation results are less 
than the symbol size. 



A. Thermal expansion 

The temperature dependence of the crystal volume at 
atmospheric pressure derived from the QHA is compared 
to PIMD results in Fig. HI Both quantum and classical 
limits are shown for the three ice phases. The QHA re- 
produces the simulation results rather accurately up to 
temperatures of 100 K. At higher temperatures anhar- 
monic effects not included in the QHA cause deviations 
of different sign in the ice Ih and III curves, while for 
ice II the agreement with the PIMD data is good up to 
the highest studied temperature. The QHA estimation 
of the zero point expansion (i.e. the volume increase with 
respect to the classical limit at T=0) of the q-TIP4P/F 
model amounts to 4% for ice Ih and II, and to 3.6% for 
ice III. 

The comparison of the QHA result to available exper- 
imental dat a 21 ' 64 at atmospheric pressure is presented in 
Fig. [5] for deuterated ices Ih and II. The overall agree- 
ment is good, in particular for ice II. At low temperatures 
there appears a negative thermal expansion in both ex- 
perimental and q-TIP4P/F results of ice Ih. However, for 
ice II this effect is not appreciable at the scale of the fig- 
ure. This behavior of ice Ih and ice II can be rationalized 
by the differences in the average Griineisen constant of 
the N translational modes of lowest frequency (see Tab. 
|TT]for values calculated for the H2O ices). The negative 
thermal expansion of ice Ih is a stringent test of the wa- 
ter model. An interesting comparison of several effective 



7 



Table II: Average of the 9iV harmonic modes obtained for the reference cells of the ice phases with the q-TIP4P/F model into 
groups of N modes. The average wavenumbers and Griineisen constants are a fingerprint of the vibrational structure of each 
ice phase. 
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Figure 6: Equation of state of ice Ih, II, and III at 100 
K. QHA results (lines) are compared to PIMD simulations 
(symbols) . All results were derived with the q-TIP4P /F water 

Figure 5: Temperature dependence of the volume of D2O ice model. Error bars of the simulation results are less than the 

Ih and II at atmospheric pressure. The QHA results (lines) symbol size. 

are compared to the experimental data (circles) of ice Ih (Ref. 

[U) and II (Ref. |H). 



potentials reveals that the negative thermal expansivity 
of ice Ih is succesfully reproduced by the TIP4P model. 
However this effect is absent in the TIP5P potential and 
only slightly visible with the ST2 model. 65 

B. State equation at 100 K 

The pressure dependence of the volume of the ice 
phases is presented in Fig. [5] at temperature of 100 K. 



The QHA shows good agreement with PIMD results. The 
sudden drop of the PIMD data of ice Ih at high pressures 
is due to the proximity of the spinodal pressure at about 
1.2 GPa where ice Ih becomes mechanically unstable and 
collapses into a high-density amorphous (HDA) ice in 
simulations with the q-TIP4P/F models Interestingly 
this limit of mechanical stability of ice Ih is also captured 
by the QHA. By increasing the hydrostatic pressure the 
QHA predicts a compressed ice Ih volume characterized 
by the appearance of soft phonons that reach a vanish- 
ing vibrational frequency at a pressure slightly below 1.4 



32 



3 

O 
| 



o 
> 



31- 



25 



1 1 


1 I 1 I 
D,0 


1 

T = 


i 1 

145 K 


Ih 


^ — _ 


— i — 1 — i — h 





- QHA 

exp 

1 , 


— i — | — 




— 1 — 
T = 


H — 1 — 

225 K. 


II 






— - 


- . i 


1 , 1 


■ 


i . - 



25 



0.1 0.2 0.3 0.4 0.5 
pressure (GPa) 



Figure 7: Equation of state of D2O ice at 145 K (ice Ih) and 
225 K (ice II). The QHA (lines) is compared to experimental 
data (circles) of ice Ih (Ref. H) and II (Ref. H). 



GPa. 

The best overall agreement between PIMD and QHA 
results is found for H ordered ice II. In the case of ice 
III a significant deviation is found at negative pressures 
of about -0.5 GPa, but for positive pressures the QHA 
provides accurate results. 

The experimental results for D2O ice at 145 K (ice 
Ih)S& and 225 K (ice ITp- are compared to the QHA 
expectation in Fig. [7J The experimental data show the 
pressure window where each phase is stable at the studied 
temperature. Note the large volume difference between 
the two ice phases. The QHA predicts reasonable results 
for both phases specially in the low pressure region of 
each phase. The main difference to the experimental data 
is seen by the slope of the V(P) curve, i.e., the QHA using 
the q-TIP4P/F model overestimates the hardness of both 
ices Ih and II (it predicts a larger bulk modulus). 



C. Bulk modulus 

The temperature dependence of the bulk modulus of 
the three studied ice phases is presented in Fig. [5J The 
bulk modulus in the MD simulations is calculated from 
the fluctuation formula of the cell volume in the NPT 
ensemble^ At any given temperature the order of in- 
creasing bulk modulus (less compressibility) is: ice III < 
Ih < II. The QHA of ice II provides results in excellent 
agreement with quantum PIMD and classical MD simula- 
tions. Also for ice Ih we find that the QHA provides good 
results. A worse agreement is found for ice III where the 
QHA predicts a bulk modulus somewhat larger than the 
PIMD results. Classical and quantum results above 100 
K are not very different from each other specially for ice 
III and Ih. The statistical error of the bulk modulus in 
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Figure 8: Temperature dependence of the bulk modulus 
of ice Ih, II, and III derived with the q-TIP4P/F mode at 
atmospheric pressure. The QHA results in the quantum and 
classical limits (lines) are compared to PIMD and classical 
MD simulations, respectively (symbols). Typical error bars 
in the simulation results are shown at 100 K (ices II and III) 
and 150 K (ice Ih). 



the MD simulations is large. The classical MD results for 
ice III have not been represented as its statistical error 
was of the order of its difference to the PIMD data. 

The pressure dependence of the bulk modulus of ice Ih, 
II, and III calculated with the QHA is compared in Fig. 
M with the results of PIMD simulations at 100 K. The 
increase with pressure of the bulk modulus predicted by 
the QHA agrees rather accurately with the PIMD data 
in the case of ice II and ice Ih. For ice III we find again 
a worse agreement between QHA and MD simulation re- 
sults. Nevertheless, the largest deviations between both 
sets of data are found either at negative pressures or at 
positive ones larger than 0.35 GPa. This value is the 
coexistence pressure experimentally found at 100 K be- 
tween ice III and the higher pressure phase ice VM- 



D. Kinetic energy 

The internal energy is given by the sum of potential 
and kinetic energy contributions. Within the QHA both 
energy terms must be identical as a consequence of the 
virial theorem. A comparison of the QHA and PIMD es- 
timations of the kinetic energy, K, at atmospheric pres- 
sure is presented in Fig. [TO] as a function of temperature. 
For ice Ih we compare also the quantum results to the 
classical limit (broken line), that amounts to ksT/2 per 
degree of freedom (equipartition principle). The QHA 
estimation of the zero point kinetic energy of ice Ih is 
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Figure 9: Pressure dependence of the bulk modulus of ice 
Ih, II, and III calculated with the q-TIP4P/F model at 100 
K. QHA results (lines) are compared to PIMD simulations 
(symbols) . The estimated error bars of the simulation results 
are shown at P = 0. 



34.4 kJ/mol. Nearly the same value is obtained for ice 
III, while for ice II we get a reduced zero-point energy 
energy of 34.0 kJ/mol. This difference implies that ice 
II is stabilized with respect to ice Ih and ice III due to 
its lower zero point energy. The comparison of QHA and 
PIMD data shows that the QHA overestimates the value 
of K by about 0.8 kJ/mol in the three ice phases and 
that this shift is nearly temperature independent. 

The pressure dependence of the kinetic energy at 100 
K is displayed in Fig. [11] Here we also find that QHA 
results are shifted with respect to the PIMD data. The 
shift remains nearly constant with the pressure, except 
for the case of ice Ih at the highest simulated pressures, 
where larger deviations are originated from the proximity 
to the limits of the mechanical stability of this ice phase. 



E. Enthalpy 

At atmospheric pressure the enthalpy of ice is nearly 
identical to its internal energy as the term PV results to 
be vanishingly small. The enthalpy, H, of the three ice 
phases at atmospheric pressure is presented in Fig. [T21as 
a function of temperature. For ice Ih we have compared 
both the classical and quantum limits of the QHA to the 
corresponding simulation results. At the scale of the fig- 
ure the agreement is very good. The QHA zero-point 
energy of ice Ih is 68.9 kJ/mol. The QHA systemati- 
cally overestimates the enthalpy of the three ice phases 
at low temperatures. In the case of ice Ih we find that 
at 50 K the QHA result is 0.7 kJ/mol larger than the 
PIMD result, while for ice II and III we obtain a larger 
deviation of about 0.9 kJ/mol. Both QHA and PIMD 
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Figure 10: Kinetic energy of ice Ih, II, and III at atmospheric 
pressure calculated with the q-TIP4P/F water model. The 
QHA results (full lines) are compared to PIMD (symbols). 
For ice Ih both quantum (q) and classical (cla) limits are 
shown. Error bars of the simulation results are smaller than 
the symbol size. 



results show that at 50 K the enthalpy of the ice phases 
increases in the order Ih<II<III. With respect to the en- 
thalpy of ice II, the relative values found for ice Ih and 
III at 50 K are -0.3 and 0.5 kJ/mol (PIMD) and -0.4 and 
0.5 kJ/mol (QHA), respectively. Interestingly at higher 
temperature the agreement between the PIMD and QHA 
data becomes better as a consequence of an error com- 
pensation between kinetic and potential energy terms. 
Thus the error of the QHA enthalpy estimation is lower 
than that found for the kinetic energy in the previous 
Subsec. llVDl 

The pressure dependence of the enthalpy at 100 K is 
represented in Fig. (T3] The term PV plays here an im- 
portant role. We have already seen in Fig. [12]that at 100 
K and atmospheric pressure the QHA overestimates the 
reference enthalpy of PIMD simulations. This behavior 
results nearly independent of the external pressure. In 
Fig. [13] the QHA result for the three ice phases lies sys- 
tematically slightly above the corresponding simulation 
results, but the overall agreement can be considered as 
satisfactory in the whole studied pressure range. 



F. Heat capacity 

The heat capacity at constant pressure is defined as 
'dH 



° p ~ I dT 



(10) 




Figure 11: Kinetic energy of the studied ice phases as a 
function of the pressure at 100 K. QHA results (lines) are 
compared to PIMD simulations (symbols). Error bars of the 
simulation results are in the order of the the symbol size. 



Figure 13: Enthalpy of the studied ice phases as a function 
of the pressure at 100 K. QHA results (lines) are compared 
to PIMD simulations (symbols). Error bars of the simulation 
results are smaller than the symbol size. 
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Figure 12: Enthalpy of ice Ih, II, and III at atmospheric 
pressure calculated with the q-TIP4P/F water model. The 
QHA results (lines) are compared to PIMD and MD simula- 
tions (symbols) . For ice Ih both quantum and classical limits 
are shown. Error bars of the simulation results are less than 
the symbol size. 



The PIMD estimation of C p of ice Ih at atmospheric pres- 
sure has been obtained from a numerical fit of the H(T) 
curve shown in Fig. 1121 The temperature derivative was 
calculated from a 4th degree polynomial fit of H(T) in 
the interval 80-190 K. The result is plotted as filled cir- 
cles in Fig. 1141 For comparison the experimental C p data 



of ice Ih of Refs. I6al69l have been plotted as open circles 
and squares, respectively. We find good agreement be- 
tween the PIMD results and the experimental data. The 
figure also includes lines showing the C p values derived 
from the QHA of ice Ih, II, and III, with an inset that 
highlights their low temperature behavior. Note that the 
QHA underestimates the PIMD data of ice Ih at temper- 
atures above 100 K. However, the QHA result appears 
close to the simulation at 80 K. At low temperatures 
(below 40 K, see the inset of Fig. [HJ the QHA for ice Ih 
shows an excellent agreement to the experimental data. 
We expect that such a good agreement will occur also in 
PIMD simulations. However, the computational cost of 
these PIMD simulations increases as 1/T, along with the 
increase in the number of beads, so that a reliable sim- 
ulation of C p below 40 K becomes computationally too 
expensive. Therefore the QHA seems to be a practical 
alternative to PIMD simulations in the study of the low 
temperature heat capacity of ices. 

The heat capacities of the three ice phases are signif- 
icantly different in the studied temperature range. In 
particular, below 50 K ice II displays the lowest heat 
capacity at the studied pressure, while the opposite be- 
havior is observed above this temperature. 

It has been previously stressed that the pressure de- 
pendence of C p in ice Ih has not been studied experi- 
mentally. Thus the empirical equation of state of Feistel 
and Wagner has been applied to extrapolate the temper- 
ature dependence of C p at pressures up to 0.2 GPa^S 
In particular the temperature dependence of the relative 
heat capacity 

AC P = C P {P) - C p (P ref = 0) (11) 
derived at P=0.2 GPa with this empirical equation of 




Figure 14: Heat capacity of ice Ih, II, and III at atmospheric 
pressure. Lines are results derived by the QHA. Filled circles 
correspond to PIMD simulations of ice Ih. Experimental data 
for ice Ih are shown by open circles (Ref. 169 ) and open squares 
(Ref. 68j). The inset emphasizes the low temperature limit. 



Figure 15: Relative heat capacity, C P (P) ~ C v (P re f = 0), 
of ice Ih as a function of temperature. The QHA results are 
shown by lines labeled by the calculation pressure (in GPa). 
The curve with symbols was derived by the empirical state 
equation of Ref. 70; at a pressure of 0.2 GPa. 



state is represented by full symbols in Fig. [TS] This ex- 
trapolation of AC P is characterized by having a rather 
small slope at low and at high temperatures. The same 
qualitative behavior is found for lower pressures. 70 On the 
contrary, the QHA predicts a different overall behavior, 
with a maximum of AC P at low temperatures (~ 25 K) 
and a conspicuous negative slope at high temperatures. 
We believe that the QHA provides a more realistic esti- 
mation of AC p than that obtained by extrapolation of 
the empirical equation of state. Our main argument is 
that the QHA represents a realistic physical model of the 
ice phase. 



V. CONCLUSIONS 

In this work we have undertaken a systematic study 
of the capability of the QHA to reproduce several ther- 
modynamic properties of ice Ih, II, and III as a function 
of both pressure and temperature. The studied pressure 
range goes from -1 to 2 GPa, while the temperatures were 
studied up to 300 K. Thus the region in the P — T plane 
where the QHA has been checked is much larger than 
the area where the studied ice phases are found to be 
thermodynamically stable. Therefore an important con- 
sideration of the present study is its generality, in the 
sense that it is not limited to a small number of state 
points. 

The validity of the QHA is restricted by the presence 
of anharmonic effects not included in the approximation. 
Thus a direct check of the QHA is the comparison to nu- 
merical simulations that fully consider the anharmonicity 



of the interatomic interactions. We have employed the 
empirical q-TIP4P/F model for this comparison. Our 
expectation that the validity of the QHA is largely in- 
dependent of the employed water model is based on the 
assumption that the anharmonicity of the q-TIP4P/F 
potential is reasonably realistic in view of the large body 
of experimental data reproduced by this model ^ 30 ' 32 

The comparison of the QHA to the PIMD simulations 
shows a remarkable overall agreement for the three ice 
phases. The best agreement has been found for ice II, 
which is an ordered phase of ice. Both crystal volume and 
enthalpy have been shown to be reasonably reproduced 
by the QHA. This fact let us expect that the QHA could 
be appropriate to study phase coexistence of ice phases 
as the slopes of phase coexistence lines are a function of 
the differences in volume and enthalpy of the correspond- 
ing ice phases (Clausius-Clapeyron relation). The QHA 
is also particularly appropriate to study the low temper- 
ature limit of certain thermodynamic properties as the 
heat capacity, bulk modulus and internal energies. The 
computational cost of PI simulations increases as 1/T, 
so that the simulation of low temperature limits can be 
prohibitively expensive by this approach. Further advan- 
tages offered by the QHA are the lack of statistical errors 
and the easier checking and correction of finite size er- 
rors. These advantages apply not only compared to PI 
simulations, but even compared to classical MD. 

The experimental heat capacity of ice Ih at low tem- 
peratures (T <40K) is reproduced with remarkable accu- 
racy with the QHA at atmospheric pressure. Our results 
lead us to expect that this good agreement should not 
deteriorate at higher pressures. We have shown that the 
temperature dependence of heat capacities at finite pres- 
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sures predicted by the QHA differs in important aspects 
from a previous estimation based on the extrapolation 
of the Feistel- Wagner equation of state for ice Hi.— Al- 
though we have not found experimental data of the heat 
capacity of ice Ih above the normal pressure, we expect 
that the temperature dependence predicted by the QHA 
is physically sound because it is based on a reasonable 
physical model of ice. 

An interesting aspect of the QHA is that it is sen- 
sible enough to predict different anharmonic behaviors 
as a function of the employed water model. The ther- 
mal expansion of ice Ih at low temperatures is predicted 
by the QHA to be negative for the q-TIP4P/F and the 
TIP4P potentials but positive (or slightly negative) for 
the TIP5P and ST2 models^ Moreover, the isotope ef- 
fect in the crystal volume of ice Ih is predicted by the 
QHA to be anomalous (as in the experiment) with a 
DFT functional, but normal with the q-TIP4P/F model. 
We stress that these differences in the QHA predictions 
are in agreement to the results of available computer 
simulations^ 



At last we mention several lines where the QHA is 
likely to provide useful results in relation to ice phases 
investigations: i) the dependence of thermodynamic vari- 
ables with system size; ii) low temperature studies of ice 
phases where the computational cost of PI simulations 
becomes increasingly large; Hi) isotope and quantum ef- 
fects in the phase diagram of ice; iv) dependence of ther- 
modynamic variables on hydrogen disorder; v) check of 
improved water models. 
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